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Symbols: 

C = specific heat of the fluid 

d = vertical distance of the tube from the heated wall 

(Figure 1 (a ) ) 

F = configuration factor 

h = convective heat transfer coefficient at the solid- 
fluid interface 

k = thermal conductivity 

L = horizontal distance of the centre of the tube from 
the left end of the heated wall (Figure l(a)) 

M = total number of grid points 

Nu = Nusselt number 

Pr = Prandtl number 

q = heat flux 

r = radial coordinate 

Zir = step size in r-direction 

Re = Reynolds number 

T = temperature 

Th = thickness of the tube wall 

u = axial velocity of the fluid 

W = width of the heated wall (Figure l(a)) 

X = axial distance 

AX = step size in x-direction 

Greek letters: 

a - thermal diffusivity of the fluid 

e = emissivity of the tube material 

u = coefficient of viscosity of the fluid 
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V = coefficient of kinematic viscosity of the fluid 

P = density of the fluid 

o = Stef an-Boltzmann cons'cant 

0 = circumferential coordinate 

A0 = step size in jJ-direction 

Subscripts: 

fb = fluid bulk 

e = at the entrance of the tube 

f = fluid 

1 = index notation for a node in r-direction 

j = index notation for a node in 0-direction 

k = index notation for a grid in x-direction 

P = value at previous x-station 

p = at constant pressure 

r = along r-direction, entirely in the fluid and wall 

of the tube 

rf = along r-direction, only in fluid region 

S = surrounding medium 

s = solid 

W = heated wall 

x = local value 

Superscripts: 

- = average value 

# r= non-dimensional value 

n = value at the previous iteration. 
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ABSTRACT 


A simple but computationally efficient finite difference 
scheme has been used to solve the problem of conjugate heat tran- 
sfer in a laminar flow of high Prandtl number fluid in a circular 
tube whose periphery is exposed to non-uniform radiation heat 
flux from a large heated wall. The present model is a simulation 
of the effect of hot refractory wall close to a tube in a pipe 
still. 


The results indicate that in a thin tube with low conduc- 
tivity, the temperature of the tube is very high which may result 
in hot spots on the tube surface* If the tube material is thick 
and has a sufficiently high conductivity, the peripheral conduc- 
tion in the tube wall alleviates this problem. Also, the temper- 
ature of the tube goes up when the tube is placed closer to the 
heated wall. 

For the wider wall the mean temperature of the inner 
wall of the tube is more as the tube receives more radiation. 
Also, it is observed that mean tube inner wall temperature incr- 
eases with the axial distance and is lower for higher Prandtl 
number fluid. 

It is also seen that thermal entry lengths are higher 
for higher Prandtl number fluid and for higher Reynolds number 
of the flow. 

The fluid bulk temperature increases with the axial 
distance and for higher Prandtl number of the fluid and lower 
Reynolds number of the flow, the fluid bulk temperatures are 
less. 


The accuracy of the present numerical solution has been 
checked by comparing it with the analytical solution of the case 
of fully developed constant heat flux condition on the periphery 
of a very thin tube. The result is found to be very satisfactory. 



CHAPTER 1 


INTRODUCTION 


1 .1 General 

This investigation presents the case of heating of a 
fluid flowing through a long circular tube which is exposed to 
a non-uniform radiation heat flux on its periphery by the proxi- 
mity of a very large heated wall situated parallely above the 
axis of the tube. The temperature of the fluid is uniform at 
the inlet of the tube where the heat transfer begins, but the 
velocity profile is already fully developed and invariant. 

This, however, limits the applicability of the present analysis 
to fluids whose Prandtl numbers are very high relative to 1 , 
for example, oils. The present analysis is for constant 
property fluids but variable property fluids can also be inves- 
tigated with little modification in the solution procedure. 

The present model is a simulation of the effect of hot 
refractory wall close to a tube in a pipe still. Frequently, 
in the case of uneven distribution of heat flux around the 
tube, hot spots on the tube surface may occur. If the tube 
material is thick and has a sufficiently high conductivity, 
peripheral conduction in the wall may alleviate . this problem, 
but for thin walled tubes, the difficulty can be accute. 

The temperature distribution in the wall is a strong 
function of the amount of radiation heat flux received at the 
outer surface of the tube and the heat transfer at the interface 
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between the inner surface of the tube and the fluid. Thus, 
the conjugate heat transfer at the interface plays an impor- 
tant role in this problem. 

1.2 Objective 

The objectives of the present investigation are (a) to 
find the temperature distribution in the wall of the tube for 
various thickness and thermal conductivity of the tube and for 
various positions of the tube with respect to the heated wall, 
(b) to get a qualitative indication of the thermal entry lengths 
for fluids of different Prandtl numbers and for various Reynolds 
number of the flow, and (c) to show circumferential variation 
of solid-fluid interface temperatures for different widths of 
the heated wall. 

1 .3 Literature Survey 

The problem of laminar flow in a circular tube with 
fully developed constant heat rate conditions axially, but with 
an arbitrary peripheral variation of heat flux, has been treated 
by Reynolds [1, 2]. In [1 ], a solution is obtained for a pulse 
of heat transferred across an elemental section of wall, and 
superposition is employed to obtain solutions for any prescribed 
heat flux. In [2], the prescribed heat flux is expressed in 
terms of a Fourier expansion and also the flow is turbulent. 

An interesting conclusion of [2] is that the effects of circum- 
ferential heat flux variation in turbulent flow are sometimes 
more pronounced than in laminar flow. The coupled effect of 
wall and fluid conduction in laminar pipe flow has been studied 
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by Faghri and Sparrow [3] and Zariffch et al [4] using finite 
difference and by Campo and Rangel [b] analytically. Davis 
and Venkatesh [6] also studied heat transfer to Poiseuille 
flow in a pipe with heat conducting section, as an example of 
their new general formulation of the solution of conjugate 
boundary value problem. Computations, however, were not per- 
formed. Luikov et al [7] presented a closed form analytical 
solution for the conjugate heat transfer problem in case of 
fluid flow in circular tube. S.V. Patankar [8] suggested the 
use of harmonic mean of the transport properties on both sides 
of an interface to treat the points of discontinuity. Therm- 
ally fully developed flow in a square duct with finite wall 
thickness was presented as an example. The method is quite 
general and can also calculate fluid flow and heat transfer in 
arbitrary geometries. Barozzi and Pagiliarini [9] have used 
finite element method to solve axi-symmetric conjugate heat 
transfer in a pipe with uniform heat flux specified on its 
periphery. The method combines a finite element treatment of 
the energy equation in the tube wall and the application of the 
Duhamel’ s theorem at the interface. 

1 .4 Salient Features of the Present Investigatio n 

The present investigation differs from the previous 

works on the following counts. 

(i) the non-uniform peripheral heat flux on the outer 

surface of the tube is actually calculated, instead 
of being prescribed. 



(ii) 

(iii) 

(iv) 


(v) 


the flow is thermally undeveloped, 
the problem is non- a xi symmetric, 
a simple but computationally efficient finite 
difference scheme has been used to solve the 
problem. 

it has the capability of handling variable 
property fluid and tube material with little 
modification in the numerical method. 
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CHAPTER 2 

PROBLEM FORMULATION 
2»1 Description of the Problem 

The physical model is depicted in Figure l(a) and 
Figure l(b). The refractory wall has a width W. The pipe, 
which has inside and outside radii r^ and r 2 » is separated by 
a distance d vertically from the wall, the tube being at a 
■distance L from the left extreme point of the wall. L and d 
together determine the position of the tube with respect to 
the wall. The wall behaves as a black body radiator at temper- 
ature Ty^. The surrounding is at temperature Tg. The upper 
portion of the periphery of the tube receives radiation heat 
flux both from the wall and the surroundings while the bottom 
portion receives radiation heat flux primarily from the surr- 
oundings. So, there is non— uniform radiation heat flux on the 
periphery of the tube (Figure l(a)), but the radiation heat flux 
is uniform in the axial direction since both the wall and the 
tube are very long in that direction (Figure 1(b)). A fluid 
enters the tube at a uniform inlet temperature T^ while the 
velocity profile is fully developed. 

So, the thermal boundary layer develops along the 
axis of the tube until the temperature profile is fully deve- 
loped. The local heat flux density at any x position from the 
inner pipe wall to the fluid is, 

q, = q,(0) = h(T, - 

where T-^j^ is the fluid bulk temperature at that x-position 
and h is the heat transfer coefficient at the given x (the 
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FIG. 1(a) SCHEMATIC DIAGRAM OF THE PROBLEM AND THE 
COORDINATE SYSTEM IN (r,cf>) DIRECTION 



FIG. Kb) SCHEMATIC DIAGRAM OF THE PROBLEM AND THE 
COORDINATE SYSTEM IN (r,x ) DIRECTION 
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method of determination of their values shown later). is 
the wall temperature at same x- station. The polar coordinates 
(r, i3) are shown in Figure l(a). 

The corresponding radiant flux density to a point on 
the outer surface where the temperature I 2 ~ ^ 2 ^^^ 

qj = q2(0) = - T^) + (1 - 

The geometric view factor fraction of radi- 

ation leaving a small area on the outer pipe surface that is 
directly intercepted by the wall; and the remaining fraction, 

1 - F^yj is absorbed by the remainder of the surroundings. It is 
described in details in Appendix A, 

The temperature distribution in the wall of the tube 
and the fluid at a particular x-location are determined by sol- 
ving simultaneously the energy equation in the wall and the 
fluid and the interface and then the solution marches in the 
x-direction. 

Following assumptions are made to model the problem 
mathematically; 

(a) The heated wall behaves as a black body radiator. 

(b) The surrounding atmosphere is transparent to thesmial 
radiation. 

(c) Natural convection from the heated wall to the tube 
surface is neglected, radiation being the principal 
heat transfer mode. 

(d) The flow is steady, laminar and incompressible. 

(e) Physical properties of both the fluid and the tube 
wall are constant. 
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(f) The tube is smooth so that frictional loss is negligible. 

(g) Axial conduction in the fluid is neglected as in the 
present case RePr is greater than 100. 

(h) Axial conduction in the tube wall is neglected as the 
tube is sufficiently long. 

(i) Viscous dissipation is neglected as it is a low speed 
flow. 

2.2 Basic Equations 

For incompressible, constant property fluid, the comb- 
ined momentum and continuity equation for fully developed flow 
becomes, 

r dr dr' p dx 


du 


with the boundary conditions, 

(i) no slip condition, u = 0 at r = r^ 

(ii) the velocity is maximum at the centre i.e. ^ = 0 at r 
Solving the momentum equation subject to the above boundary 
conditions, the velocity profile obtained as follows: 

u = Average velocity (axial) 


= 0 . 


%ax _ 1 dp 2 

= ~ 2 ~ - Sp dx 


as. 


u 


max 


= u (r = 0) = - i 


1 


^ _2 


4p dx 1 


Inserting u in the expression of u(r), the velocity profile 
gets the form 



9 


u(r) 

whe re 
Energy equationst 


2u u^^(r) 


u^ir) = 1 - ( 7 “) 


( 2 . 1 ) 

( 2 . 2 ) 


(i) For solid region; 7-5 + r Ir "*■ 2 .0,2 


_ 1 3T , 1 3^ 

3 r j,2 

with boundary conditions. 


= 0 


( 2 . 3 ) 


at r = r^ 


£L 

s 3r 


_ k ^ = - = radiation heat 


flux falling on the upper surface 

s 


|1 = 1- ,5(0) = f [F2,(0 )(t: - T^) + t1 - 


( 0)1 




3 T 


. at r - r2 , 
where F2(^ = F2^(0) = 


F- Cf» Tw + (1 - 

s 


^2 =T| = 


(ii) Fluid region; 


configuration factor, determined by Hottel 
Cross-string method (see Appendix A) 

T( 0 ) 


Energy equation; u ^ 


. « - , r^T 

♦ ^ ?)X ^^3^ ^ r^ 30^ 


3T 


A O i ^ A 

with boundary conditions, at r = 0, - 

at X = 0, T = Tg 

At. the interface, the equality of temperature and heat flux 
exists, so that, 


( 2 . 4 ) 

( 2 . 4 a) 
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at r = ri , T |^ = T 

and - k ^ = - k - ^ 

s 3 r 1 f 3r 

s 


(2.5) 


2.3 Fluid Bulk Temperature 

For the calculation of the convective heat transfer 
coefficient h(x) from the wall to the fluid along the axis of 
the tube (x-direction ) , it is necessary to define a character- 
istic temperature of the fluid commonly known as 'bulk tempera- 
ture* or 'mixing cup temperature'. The reason is that in a 
tube flow there is no easily discernible free-stream condition 
as is present in the flow over a flat plate. For most tube or 
channel flow heat transfer problems the topic of central 
interest is the total energy transferred to the fluid in either 
an elemental length of the tube or over the entire length of 
the channel. At any x position, the temperature that is indi- 
cative of the total energy of the flow is an integrated mass- 
energy average temperature over the entire flow area. The 
numerator of equation (2.6) below, represents the total energy 
flow through the tube and the denominator of the same equation 
represents the product of mass flow and specific heat integ- 
rated over the flow area. The bulk temperature is thus repre- 
sentative of the total energy of the flow at a particular 
location - (x). 

Mathematically the bulk temperature for the fluid is 


defined as, 



1 1 


Thus 


2-rt 

/ / PC uTrdrdj3 

o o P 

2n 

/ / P C u r dr d0 

o o ■ 


Tfb = 


( 2 . 6 ) 


The heat flux from the wall of the tube to the fluid at a given 
X— location and circumferential location 0, Q')(x» 0) defined 
as. 


q^(x, 0) = -[- kg 


] = k 


ir=r-, , s 


O. 

s 3r 


(2.7) 


■r=r^ ,s 

and the average heat flux at a given x is defined as, 

1 

= h -' "ii 

o 

The average wall temperature of the tube at a given x is expres- 
sed as, 


(2.7a) 


2Ti 


T,(x) = h ' h 


(2.8) 


where T^ = T^(x,, 0) = Interface temperature = (2.8a) 

Local heat transfer coefficient, h^^ can be found by the foll- 
owing heat flux balance equation; 




(x) = - “^fb^ 


or- 


hx = 


q-i ( x) 


- Tfb 


(2.9) 


The local Nusselt number, Nu is found by, 
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Nu^ 


h^(2r^) 


( 2 . 10 ) 
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CHAPTER 3 
METHOD OF SOLUTION 

3,1 Finite Difference Approximations 

The energy equations of the solid and the fluid along 
with the compatibility conditions at the interface are discre- 
tised using finite difference approximations as described below 
A typical orthogonal mesh system for using the Finite 

Difference approximation is shown in Figure 2. 

The following finite difference approximations are used 

3^T 3_T 

(i) Using central difference formulae for the terms ^^2 » 3r 

and at the point (i, j), 

30 ^ 


3^T 

1 

^”“'i+1J 

"■ ^"^1,3 ”^1-1 » j 

I? " 

(^r)^ 


1 


- 2T^^^ + 

30^ 

(A0)2 

3T 

1 

*^"^1+1 , j 

“ Vl,j ^ 


2C^r) 


where Ar = uniform step size in r-direction 
A0 = uniform step size in 0-direction 
i = index in r-direction 
j = index in 0-direciion« 

(ii) Using backward difference formula for the convection term 

^ at a x-station, k+1 , 

3X 

3T , i ,k-H - i.j tk 

§7 = 
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where T. . = Temperature of the point (i, j) at the previous 

X- station. 

Also term can be discretised by backward difference formula 
3r , 


as, 




Ar 


3.1.1 Discretization of the E nergy Equation for Solid 

The equation (2.3) subjected to the boundary condition 
(2.3a) can be discretized as follows, 

4.. + i ^ = 0 

Equation (2.3); + r 3r ^ ^ 

Boundary condition (2.3a); at r = 12 ; 

I? = F ■''w ■ ''2V(^'‘'s - '^2 ^ 

s 

Refer to Figure 2. 

Let, Nj. = Number of grid points in r-direction for the entire 
fluid-solid region 

N,, = Number of grid points in 0-direction, both for solid 

0 

and fluid region 

N = Number of grid points in r-direction, only in the 
rf 

fluid region. 

As x-direction is taken Infinitely long, no specific number of 
grid points are taken for this direction. By marching procedure 
the solution can be continued upto any length and it can be 

terminated as per requirement. 

The uniform step sizes are. 
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Ar, 


(r2 - 


A0 = 


2n 


N 


0 


and 


Ar^ = 


TN 


rf 


'1 

^=rTJ 


Discretizing equation (2.3) by finite difference approximation, 
as 


( A r 


^2CVl,j “ ^^i,j ^i+1,3^ 2A-rg '^’^i+1,j " "^i-IJ^ 


^ ^ [T. . 1 - 2T. . + T. . J = 0 

1»J+' 


Ar 2 


(A0) 

Simplifying the above equation, 

(21^ - 

Ar 2 , 

+ (Zif + Ar3)T.^,_j = - 2(^) [Tij., + 

Writing the above equation in the standard tridiagonal matrix 
form. 




(3.1 ) 


A. = 2r^ - r. Ar^ 

1 1 IS 

o Ar 2 

Bi = -[4r^ + 4(-2^) ] 

Ci = 2r? + r. Ar 

^ JL w 

Ar 2 

Di = - 2(s^) + '^i,j+1^ 


where 
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The above equation is valid for any interior point (i, j) of 
the tube wall except the boundary points, i.e® for all nodes 
numbering from i=N^^ + 1 toi=N^-1 and in 0-direction 
all points except the points at 0 = 0 and 0 = (2Tt - 0 ).i.e. 

j = 1 and j = 

More correctly saying, the equation (3.1) is valid also 
for j = 1 and j = N0, but for computational purpose the equa- 
tion is slightly modified for all points with j = 1 and j = N0, 
as described below. 

For j = 1, the point (j-1 ) physically represents the 
point '1^0' • Therefore the equation (3.1 ) is modified by repla- 
cing (j-1 ) by N0 for j = 1 and ( j+1 ) by 1 for j = N0, both 

changes occur in the expression of 

Rewriting equation (3.1) for different j's the general 

foiTO becomes, 




(3.1) 


With = 2r^ - r^ Ar^ 


o 2 

= - [4 r^ + 4(-^-^) ] p for all j 


and 


Cl = 2ri + ri Ar^ 

Ar 2 


Ar 2 

= -* 2(-^) + ■^i,2^ for j = 1 


'0 

Ar 2 T . 

s \ Tt j. t j for j = N 


" ■“ ^^A0 ^ ^”^1,1 ’’’ ^i,N0-1 


0 
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However the equation (3.1 ) is only valid for interior points 
of the tube wall. 

To get the discretised energy equation for the upper 
surface, i.e, i = (r = the boundary condition (2.3a) is 

used to modify equation (3.1) as described below. 

Boundary condition (2.3a) is, at r = r 2 (i.e* i = 

9T ^ rp + h F ) ] 

S 


where T 2 = temperature of the tube at r = r 2 , this is repre- 
sented by T» 

Putting T in place of I2 discretlsing equation (2.3a) 

T. - T. ^ ^ r_ _4 /. ^ \ ^4 

"s 


. (1 - - it,} ’ 


but as if j term introduces non-linearity in the discretized 
equation it is necessary to linearize this term using Patankar’ s 
approach [l0] to linearization of the source term. 

According to this approach, the non-linear term should 
be linearized in such a manner so that fast convergence is 
attained. It also suggests that when the coefficient of the 
transformed linear term is equal to the slope of the original 
nonlinear term the convergence is fastest, however this is not 

true always. 

Using this approach term is linearized as. 


if . 


= 4(T^ 3(^n 








where T? . = Temperature at point (i, j) in the previous 


iteration. 
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Putting this in the discretized equation of the boundary condi- 
tion (2.3a), the equation becomes, 


T _ T. . 

^i-fl.i 1-1, j 

2Ar^ 


oe 

icT 




Simplifying it, 

Ti+I.j = ^ [''2 WjTw + (' - 

Putting this expression in equation (3.1 ) the discretized energy 
equation for i = is obtained as follows. 

Vi- 1.3 ^ Vi. = °l 

4r^ = 4r| (as i = =» = r^) 

O 2 


whe re A. = 


Bi 


c. = 0 


D. 

1 


where 


Ar. 2 T 

- S2[^2W.^W + ' ^2W. ^"''s ^ 

J J 

= ^ s^^^2 ■*■ ^2 


The above equation is valid for i = only, here also 
D. is modified for j = 1 and j = as before, i.e. replacing 
(j_1 ) by N^, for j = 1 and (j+1 ) by 1 for j = N^. 
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3.1.2 Discretization of the Energy Equation for Fluid 

The energy equation for fluid (Equation (2.4)) is 
changed by inserting the expression of u(r) into it as below. 


^ ^ 3x “ “f L - 2 r 3r 2 ^.2^ 

ox r op 


or. 


2u 


u*(r) 


_ 1 ? 1 4 . i il . ifl 

3x - 3^2 r 3r ^2 3^2 


Discretizing this equation using finite difference 
approximation for any point (i, j ) in the fluid region and at a 
X- station (k+1 ), 


2u u^ T,. . 
‘f 

1 1 


T.. 


1 1,1, k+1 i,i,k ^ -1 Ft + T 1 

— . ' * l-VlJ - + ^i,j+lJ 


k+1 


For simplification the two dimensional notation for T 

is used, replacing ^p. ^ '^i.j^k+l '^i,j* 

Simplifying the above discretized energy equation for 
fluid and writing it in the standard form of tri-diagonal matrix, 
as below, 

Vi-i,j * Vi,j = °1 


where 


% 


2ri - Ar^ r. 


^ Ar- 2 o un 
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Ci 

= 2rf + Arj 

^i 


Ar,. 2 


^i 

1 

II 

[t 


4u(Ar^)^ 



Ax 



•. . J - C u f r? Tp 
itO+1 1 1 P. 




Equation (3o2a) is valid for i = 2 to i = - 1. As 

in the case of solid region, this equation is also modified 
for j = 1 and j = N0 and the modification takes place in the 

expression of D^. 

The boundary condition (2.4a) is discretized as below. 
3T 

For i = 1, i.e. r = 0; ^ = 0 

T • - T - 

discretizing it; ~ ^ 

(3.2b) 


or, 




Combining equations (3.2a) and (3.2b) the general form of the 
discretized energy equation for fluid is obtained as below, 

(3.2) 


whe re 


^Vl,j Vi,j * Vi-H.j = °i 


r . A r , 


Ai = 2r^ - x ^ 


B. 


Ar - 2 0 * T 

= - [4r? + + C: r^ ] 


^i = 2r. +r. Ar^ 

2 


Di = - - C uj r2 T 

1 f J 


valid for i = 2 to i = ^ 
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and also for i = 1 , Aj. = 0 


B. = - 1 

Ci = 1 

D. = 0 


Equation (3.2) is modified for j =1 and j = N, 


3.1.3 Discretization of the Compatibility Conditions at the 
Fluid- So lid Interface 

The compatibility conditions (equation (2.5)) are, 


at r = n , 


's 3r 


"f I? If 


T|s = 


Energy equation for solid, equation (2.3) 


3^2 + r 3r %2 3 f ,2 - 


and energy equation for fluid, equation (2.4) 


2u 


u*(r) 


1 ai , 1 


ajr _ 4 . ^ 

3r^ ^ 9)^^ 


Refer to Figure 3. 

From Taylor's series expansion for point (i-1,j), 




3 T' 


A r? ^2- 

+ (^) 

i,j,f 3^ 


Solving i^—k) from the above equation 

3r^ i,j,f 


(i!l) 

V 2-' 


3^ ijj»f \ “ 

At the interface uf = 0, also r^ = r^ 


-r^2 1 i - A i + ^ 

)2 1-1.0 1,0 ^ i,j,f 
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FIG. 3 FINITE DIFFERENCE NODES NEAR THE 
INTERFACE 
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Putting all these in the energy equation for fluid, 
i.e. in equation (2.4), 


0 = 


[Tj. 1 1 - Tj . + Ar^cfl) ] 


1 i,j,f r^j 30 i,j,f 


1 /3^T' 


Simplifying it, 

(il) = [T _T ■ jfl) ] 

^ - 2r, Ar. + (Ar.)2^ i.j •> 




•1 “"f 


2r:j^ 30^ i,j,f 

(3.3a) 


Expressing T. . . in Taylor's series as, 

1+ i f j 2 

T - T + At 

^i+1,j “ ^i,j ^s^ar^ . **■ 2 . 

i>J»s i,j,s 


thereby, 

(A) 

3r i , j , s 


- T. , - r3(|2) ] 


(Ar.)^ 


lfJ,S 


Using this in the energy equation of solid and simplifying it, 
the (fi) 


g^y term is expressed as below, 


i»J,s 


( 31 ) 

' a-n/ 


2r 


1 


i,j,s 2r^ Ar^ - (Ar^) 


(Ar )‘^ -2_ 

5 ["^1 + 1 i ■" i " ^ ^ — 2^ ^ 

2 i+1,j 1,1 2r^ 30^i,j,s 


(3.3b) 

Ultimately two expressions are obtained for the temperature 


gradient at the interface for fluid and solid regions as, 

2r^ rr T / 3^ 

,(2r^ + Ar^) i-1,j i,j 2rJ 0^i,j,f 


(3l) — 

^ar-'. . - Ar.^ 


2 ^ 


and 
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(— ) 


i,J»s 


2r^ r T T 

Argt2r^ - Arg) ^ ^i,j 


^ (A) 

2r^ 30^ i,j,s 


] 


The compatibility conditions (equation (2.5)), i.e. Tj^ TJ^ 


at r = and 


k 21| = k. Ill at r = r, 

s 3r I f art - 1 

S I 


and (■^) will be equal and both of them have 

same expression because the temperature of fluid and solid are 


* ^2_ 

(H) 


3^T 

the same at the interface and s is expressed as, 

30 


2^ ^2_ 

30 i,j,s 30 r» 3 *r-r^ 


= iA) 


(a0)2 


Using the heat flux equality condition, i.e. 

ir / ^T \ _ If ( 3T \ 

^ i,3,f 

and expressing (|i). . and (|I). . as in equations (3.3a) 

and (3.3b), the following equation is resulted 

k Tt t + ( ^^X ) 

IT^ • / a *i+1»j ” 2r^ 30^ i,j 

*^f 2ra Ar^ - (Ar ) ^ 


■1 '"*s 

2^1 


-Ft. . - T. a 3: 

2 ^ 1 f J ' »d 


(Arf) ^ 


' / * 

2r, Arj + Zr^ 30 

Using the expression of (|^), . in the above equation and after 

necessary simplifications the^standard discretized equation for 
the interface is obtained as below. 
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■"i 

Vi 

J + Vi,: + 

where. 


= 2r2 = C. 


Bi 

- - (2r^ R + 2r^ + RR) 


Di 


where 

R 

k,, 2r^ Ar^ + (Ar^)^ 

^f * 2r,, Ar^ - (^^g)^ 

and 

RR 

R ( ^^g)^ + ( Ar^)^ 

(A0)''^ 


(for i = Nj-f) 

(3.3) 


Equation (3.3) is valid for all j’s except j = 1 and j - Nj^.* 
For these two j-stations the equation is modified by replacing 
(j-1) by for j = 1 and ( j+1 ) by 1 for j = 

3.1.4 Simplified Equations for 


2n 


From equation (2.6), - 


/ / u T r dr d0 

o o — 

2Ti 

/ / U r dr d0 

o o 


as 


PC = Constant for the fluid. 

Putting the expression of u(r) = 2u . u^^(r), where 
u^^(r) = 1 - (~)^, ^den simplifying the above expression for 

T the following expression is obtained, 
f b 

T = ^ u*{t) T r dr d0 ^3.4) 

0 0 

The heat flux, q-j(^» 
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- T. 




= k. 


r=r^ , s 


‘Nrf+1J N^fJ 

Ar. 


(3.5) 


3.2 Numerical Solution of the Discretized Equations 


The set of discretized equations to be solved are 
listed below. 

Equation (3.2), j + ®i^i,j ^i’'^i+1,j °i 


where 

2 

= - r, 

Arf 

®i 

1 1 

1 

11 

Ar 

^i 

= 2r? + r. 

Arf 

°i 

Ar 

= - 

2 


and 


for i = 2 to i = 1 


^i 


®i 




0 

- 1 
1 


, for i = 1 


= 0 


Equation (3.3), \’'^i_1,j Vi,j ^i^i+1,j " °i 


where = 2r^ = 

2 „ ^2 


B. = - [2r,R + 2r7 + RR] 

1 1 1 
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and 




R = 


RR = 


kg 2r^ Ar^ + ( Ar^)- 


2r. Ar„ 
1 s 


(Arg)- 


R( Arg)^ + ( Ar^)^ 

( A0)2 


Equations (3.1) and (3.1a) combinedly written as 


Vl-1J 


(3.1b) 


where 


= 



Bi 

r= 

o Ar 2 


Ci 

= 

2 

2rr + r. Ar 
i i s 


°i 

== 

Ar 2 

and 


= 

3r2 


h 

= 

- [ 482 ( 1 ?^^)^ + 4r| 


^i 

= 

0 


°i 

=r 

Ar 2 


for i 
to i 


- 1 


Ar 2 


^iJ+1^ 

" ^2 ^^2W.*^W + - ^2W. ^’’•S ^ 


where = 


2 ere 


^ Arg( 2 r 2 + Arg), valid for i = 


All the above equations are to be modified for j - 1 and j - N^j 
the modifications will be only in the expression of D^. For 
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j = 1 » (j - 1 ) is to be replaced by and for j = N 0 , ( j + 1 ) 
is substituted by 1 , 


Equation (3.4), T^]^(x) = 

Equation (3.5), q^ (x, 0) 


2 2 % 

— 75 / / u*(r) T r dr d0 


- T. 


Tir^ o o 


= k. 




Equation (2.7a), 
Equation (2.8), 
Equation (2.9), 


Equation (2.10), 


1 


2u 




T,(x) = 


271 

/ 1 , d 0 


271 : *1 

O 


\ = 


Nu 

X 


q-i M 
- ’■fb 


The set of discretized equations (3.2), (3.3) and (3.1b) for 
fluid, interface and solid are to be solved by a suitable nume- 


rical technique. 

In the present analysis a line-by-line method, sugges- 
ted by Patankar [l0], is used. The line-by-line method helps 
to solve the set of simultaneous linear algebraic equations by 
Tri-Diagonal Matrix Algorithm (TDMA) method. 

The line-by-line method is described in Appendix-B in 


great detail. 
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3.3 Algorithm for the Numerical Solution 

To solve the discretized equation (3.2), (3.3) and 
(3.1b) to get the temperature at a particular x-location, the 
following steps are taken. 

(a) The computations are stated at the location x = Ax 
using the known inlet temperature at the entrance of 
the tube. 

(b) The r-direction lines are selected to initiate the 
solution procedure by line-by-line method. 

(c) All the temperatures are assumed initially, both in 
the fluid and solid region. 

(d) Start solving by TDMA with the line 0=0. 

The temperature at all the grids of this line are solved. 

(e) Go to the next line along 0-direction (i.e. at 0 = A0) 
and solve the discretized equations to get temperature 
at all nodes on that line. 

(f) This procedure is followed one after another for all 
the r-direction lines. 

(g) The steps ’d' to ’f are repeated till the process 
converges, i.e. the temperature at any node does not 
change further with iterative procedure. 

After getting the converged values of temperature, 
they are stored for that particular x-station. Then the fluid 
bulk temperature, the average heat flux, the local values 
. convective heat transfer coefficient and the Nusselt number 
are calculated following the equations (3.4), (3.5) and (2.7a) 
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to (2.10)i The integrations involved are done by Simpson's 

rule of numerical integration. 

The procedure then marches in x-direction and continued 
till the fully developed situation is reached. The fully deve- 
loped situation is indicated when there is no considerable 
change in Nu^ along x. 

The accompanying Flow-chart (Fig. 4) represents the 
algorithm pictorially. The input data to the program and the 
listing of the computer program are given in Appendix-C and 
Appendix-D respectively. 
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CHAPTER 4 

RESULTS AND DISCUSSIONS 

Figure 5 shows the temperature distribution in the 
tube wall for different conductivities of the tube material at 
the circumferential location, 0=0® (highest point on the 
tube). The temperatures in the tube wall and also the temper- 
ature gradient are higher for low conductivity material than 
for high conductivity material. This is because for low conduc- 
tivity material the resistance to heat transfer is high and 
hence higher temperature-gradient. Also, because of less heat 
transfer to the fluid, the inner wall is cooled less by the 
flowing fluid which results in raising the temperature of the 
wall. This suggests that in case of low conductivity tube 
material there is a chance of burn out of the tube at the point 
of highest temperature which in turn results in leakage of the 

fluid. 

In Figure 6, the effect of wall thickness on the temp- 
erature distribution in the tube wall is depicted. It is 
observed that with the decrease of tube thickness, the temper- 
ature gradient decreases and the temperature level in the tube 
increases. So. a thin tube with low conductivity is not 
desirable as it may give rise to hot spots occurring due to 
high temperatures on the surface of the tube. The nature of 
.the curves satisfy the physics as with the reduction of thick- 
ness, the thermal resistance to heat flow decreases and hence 
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06.5 RADIAL VARIATION OF THE TEMPERATURE IN THE 
TUBE WALL FOR VARIOUS CONDUCTIVITIES OF THE 
TUBE MATERIAL AT ^ =0* 

(d = 0-1 ,Th:0 02,W = 10,Pr = 51,Re=2000 ) 




TEMPERATURE (®K ) 
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NODAL POINTS IN THE TUBE WALL 


FIG. 6 TEMPERATURE VARIATION IN THE WALL^ 

for different thickness of the tube at 

<^ = 0- (d = 0.1,W=1-0,ks = 48 , Pr=51, Re-2000) 
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decreases in temperature gradient. Also most of the radiation 
heat input is used to raise the temperature. 

Figure 7 shows the effect of vertical location (d) of 
the tube with respect to the heated wall on the temperature 
distribution in the tube wall. As expected, the temperature 
level shoots up when the tube is placed closer to the wall as 
it receives more heat from the wall. However, the temperature 
gradient remains independent of the location of the tube (d) 
as the tube thickness and thermal conductivity are same for 
all these cases. 


Figure 8 shows the circumferential inner tube wall 
temperature distribution (fluid-solid interface temperature) 
at a particular x-looation for different widths (W) of the 


heated wall when the tube is placed eciuidistant from the two 
edges of the wall. As one would expect, the curves show symme- 
tric nature about f) = 180" as the radiation heat flux between 
00 and 180" is same as that between 180" and 350". However, 
the curve for bigger width of the wall is placed above that 
for smaller widths of the wall. Also it is more flattened. 

This is quite expected, as the tube receives more radiation 
when the heated wall is more wide. The flattened nature is 
because the circumferential radiation heat flux is more uniform 
in nature in case of larger width (W) of the heated wall. 

Figure 9 indicates the plots of mean inner wall 
temperature against the axial distance (x) for different 
Prandtl number fluids, for a constant Reynolds number. It is 
observed that the mean wall temperature is lower for higher 



temperature (“K) 
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FIG. 7 


RADIAL TEMPERATURE DISTRIBUTION IN THE TUBE 
WALL FOR DIFFERENT VERTICAL PO^mON OF THE 
TUBE (d ) V/ITH RESPECT TO THE HEATED WALL 

ahfo02> = 10.ks-A8,Pr=51.^ 



SOLID FLUID IKTERFACE TEMPERATURE (“K ) 



(DEGREES ) 


MEAN WALL TEMPERATURE , T, (“K) 



AXIAL DISTANCE, X (m) 


FIG.9 


AXIAL DISTRIBUTION OF MEAN WALL TEMPERATURE FOR 
DIFFERENT PRANDTL NUMBER (Pr ) AT A CONSTANT 
REYNOLDS NUMBER 
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Prandtl number fluid and the temperature gradient near the 
entrance of the tube is more in case of low Prandtl number 
fluid. This is because for low Prandtl number fluid, the 


thermal diffusion is much more as compared with high Prandtl 
number fluids. But, as the axial distance increases i.e, at 
larger x, the temperature gradient is not affected much by 
Prandtl number as the boundary layer thickness grows along 
the x-direction. 

Figure 10 shows the local Nusselt number, Nu^^, varia- 
tion along the axial distance, x for various Prandtl number 
fluids, Reynolds number remaining same. As expected, Nusselt 
number is very high at the entrance of the tube and then it 
decreases monotonically along x until it becomes almost cons- 
tant far downstream. Local Nusselt numbers are less for 
lower Prandtl number fluids. This is because lower Prandtl 


number fluid have higher thermal diffusivity and hence convec- 
tive heat transfer is less. So, the heat transfer coeffloient 
is also less which in turn implies that Nusselt number is less. 
Also, it is seen that thermal entry lengths for higher Prandtl 
number fluids are much more as they have lower thermal diffu- 
sivity which means that thermal boundary layer develops slower 


Figure 11 shows how local Nusselt number, Nu^ changes 
along 'x' for various Reynolds number, Prandtl number rema- 
ining same. For higher Reynolds number of the flow, Nusselt 
numbers are larger. This is because, higher Reynolds number 
means higher average velocity of the flow and hence more heat 
is carried away by the fluid and since Nusselt number is the 



LOCAL NUSSELT NUMBER (Nu 
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FIG. 10 AXIAL DISTRIBUTION OF LOCAL NUSSELT NUMBER FOR VARIOUS 
PRANDTL NUMBER (Pr) AT CONSTANT Re 



LOCAL NUSSELT NUMBER (Nux) 


2 



FIG.11 


AXIAL DISTRIBUTION 
REYNOLDS NUMBER 


OF THE LOCAL NUSSELT NUMBER FOR DIFFERENT 
(Re ) OF THE FLOW, AT A FIXED VALUE OF Pr 
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Figure 12 shows plots of fluid bulk temperature 
along the axial distance (x) for various Reynolds numbers, 
Prandtl number remaining same. The plots indicates that 
fluid bulk temperature increases with the axial distante, x. 
For higher Reynolds number of the flow, the fluid bulk temper- 
atures are less than that for lower Reynolds number. This is 
because for higher Reynolds number, the average velocity of 
the flow is higher which results in more convective heat trans 
fer giving rise to lower bulk temperature of the fluid at a 
given x-location. 

Figure 13 shows plots of fluid bulk temperature 
along the axial distance (x) for various Prandtl number 
fluids, Reynolds number of flow remaining same.- The plots 
indicate that fluid bulk temperature increases with axial 
distance, X. Also, for higher Prandtl number, fluid bulk 
temperatures are lower. This is because for higher Prandtl 
number fluids thermal diffusivity is lower which in turn 
results in higher convective heat transfer and hence lower 
fluid bulk temperature at a given x-location. 

Finally, the accuracy of the present numerical solu- 
tion is checked by comparing the Nusselt number for the fully 
cj 0 '\^glopgd condition when constant heat flux condition is 
imposed on the periphery of the tube and conduction in the 
tube wall is neglected, i.e. the tube is considered very very 
thin. Nusselt number at a distance of 2000 meters along the 
tube is found to be 4.5 which compares very well with 4.364, 
the Nusselt number for fully developed constant heat flux 
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FIG 12 AXIAL DISTRIBUTION OF FLUID BULK TEMPERA- 
TURE FOR VARIOUS VALUES OF Re OF THE FLOW 




45 



FIG. 13 AXIAL DISTRIBUTION OF THE FLUID BULK TEMPERA- 
TURE FOR DIFFERENT VALUES OF PRANDTAL NUMBER 
( Pr ) 
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condition as obtained by analytical solution. Since, theore- 
tically thermally fully developed condition occurs at infinity, 
it is not numerically possible to go to that extent because 
of limits on computation time. As the local Nusselt number 
for uniform heat flux condition obtained numerically is very 
much closer to the value obtained by analytical solution, the 
present solution for non-uniform heat flux condition can be 
called accurate enough. 
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CHAPTER 5 
CONCLUSIONS 

A simple but computationally efficient numerical 
scheme has been used to solve the problem of conjugate heat 
transfer in a laminar flow of high Prandtl number fluid in a 
circular tube whose periphery is exposed to non-uniform 
radiation heat flux from a large heated wall. The numerical 
results obtained agree very well with the physical expecta- 
tions. Also, the accuracy of the present numerical solution 
has been checked by comparing the Nusselt number obtained 
numerically for thermally fully developed condition with that 
obtained analytically for the idealised case of constant heat 
flux condition on the periphery of a very thin tube. Also, 
this analysis can be modified to take into account variable 
property fluid and the tube wall. The present analysis is 
for infinite length of the tube, but it can also be extended 
to the tube of finite length with necessary change in boundary 


conditions. 
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appendix a 

DETERMINATION OF THE CONFIGURATION FACTOR 

A. 1 Definition 

The configuration factor is also called the geometric 
view factor. For any two given surfaces, the orientation 
between them affects the fraction of the radiation energy 
leaving one surface that strikes the other surface directly. 
Therefore, the orientation of the surfaces plays an important 
role in radiation heat exchange. The physical significance of 
the view factor between two surfaces is that it represents the 
fraction of radiative energy leaving one surface that strikes 
the other surface directly. 

A. 2 Hottel*s Cross-String Method Jill 

Consider an enclosure, as shown in Figure A-1 , consis- 
ting of four surfaces which are very long in the direction 
perpendicular to the plane of the figure. The surfaces can 
be flat, convex or concave. To find the view factor, F-|_2» 
imaginary strings are assumed, shown by dashed line in the 
figure, they are tightly stretched among the four corners A, 

B, C, D of the enclosure. If ^ d = 1 ,2^3,4, 5,6 )denote the 
lengths of the strings joining the corners, as illustrated in 
Figure A- 1 . According to Hottel, the view factor F ^_2 
expressed as. 
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FIG. A-1 DETERMINATION OF CONFIGURATION FACTOR 
BETWEEN THE SURFACES A, AND A2 BY 
HOTTEL'S CROSS-STRING METHOD 



F. _ = crossed l ength - Total uncrossed length 

2 X Length of surface '1' 

(1.5 + L^) - (L^ + L^) 

~ 2 l Length of surface M’ 

1 approximated to L^ 

"^^he Equations for Determining Configuration 

For the present case the configuration factor is 
found between any point of the outer surface of the tube and 
the large heated wall placed over the tube. Each point on the 
upper surface of the tube receives radiation heat flux from 
the entire plate. 

Figure A-2 is considered, the tube is placed at a 
distance d vertically below the large wall. Both the wall 
and the tube are infinitely long in the direction perpendicular 
to the plane of the paper. 

In Figure A-2, A is point of tangent on the outer 
circle from the right end point N, of the wall. Similarly A' 
is the touching point of the tangent drawn on the outer circle 
from the left end point, N' , of the wall. 

Clearly the region below the points A and A’ do not 
get any radiation heat flux directly from the hot plate. 

Thereby the configuration factor in this region is obviously 
zero. So, the location of points A and A' is most important, 
however it is to be noted that when the tube is placed verti- 
cally below any edge then point A or A' will be shifted to 
point C or C respectively. Therefore in the Figure A-2 
angle a cannot be less than P and limiting value of a is P, 
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however if the top point of the tube, P is placed not directly 
below the plate then a can be less than p. 

To find out P the triangle MCO is considered, 

MO = d + r^ 

OC = T 2 == Outer radius 


MC 



as MCO is a right angle triangle and 
I MCO = 90° 


or MC = Ad + - r^ = '^d^ + 2d r 2 


As, Z.MOC = P , . tan P = 

. '^d^ + 2d r^ 
p = tan“^ (• j ) 


OC 


+ 


2d r. 


(A.1 ) 


Consider the triangle OAB, 

/.AOB = ( 180 ° -a) 

/.OAB = 900 

= 1800 _ [/aOB + /.OAB] = a - 90° 

In the right angle triangle OAB 

^ = sin 0/ = sin(a - 90o) = - cos a 

OB 1 

OB . = - OA/cosa = - r 2 /cosa 

.*. MB = MO + QB = d + r2 - ^2 
MN = W - L 



54 


MN 4. ^ 

* ’ MB ~ = tan(a -90°) = » cota 

W - L 

’ a + - rj sec cc = " 

or, (W - L)sina = _ (d + 42)0050 + 

Squaring both sides of the above expression and simplifying it 
the following expression is obtained 
2 

A cos a - B coset + C = 0 , 

where A = (W - L)^ + (d + ro)^ 


B = 2r2(d + r^) 


C = r2 - (W - L)‘ 


B + - 4AC 


cos a = 


2A 


B + D 
2A 


D = - 4AC 


sina = '^1 - cos 


= ^Aa^-Cb + d)^ 


when, cosa = 


B + D 
2 A 


= . (b - D)'-^ 

B - D 

when, cosa = 


or. 


> /4A^ , (B - D)^ 
B"^rD 

W .. (6 4 Djl 
B + D- 


tana 


(A.2) 
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Similarly, to find out 6, the triangle OA'B' is 
In triangle OA' B' , 

^OB>A' = , Zb'OA’ = (1800 - 6) 

. . ©2 = 1800 _ [900 + (1800 - 6 )] = (6 - 

OB' = OA' /sin = _ r2/cos6 

From triangle, MB'N', 

MB' = MO + OB' = d + r2 - r2 sec6 
MN' = L 


tan ©, 


MN' 

MB' 


(d + r2 - ^2 ^ ^ 


or, - cot 6 = X— j- 

* d + r2 - r2 sec6 

Solving this equation for ’6’, 

Aa^ - (B^ + 


tan 6 = 


B, + D, 


or. 


- (B^ - D^)- 


Bi - D, 


where, 

A 


B 


1 


1 




= /(B^ _ 4A^c”) 


= L 4- (d + Tr) ) 


2r2(d + r2) 




Y = 


considered. 


900) 


(A.3) 


Also, 


2it - S 


(A.4) 
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From the equations (A.1), (A.2), (A.3), (A.4) the values of p . 
a, 6 and y are obtained, which play important role for the 
calculation of configuration factor. 

Next step is to express the configuration factor, 

F^w as a function of '0'. Refer to Figure A-3. 

Let, the configuration factor for the elemental arc XY 
on the upper surface is required to be found. The length XY 
is sufficiently small such that in the limiting case this 
element turns to a point and the expression for F 2 ^y this 
element becomes the value of F^^ at the point P. 

Points X and Y are joined with points M and N by four 
strings, two crossed, two uncrossed, as illustrated in the 
Figure A-3. Here, XN and YN* are the uncrossed strings and 
XN' and YN are the crossed strings. 

By Hottel’s Cross-String method, as described earlier, 

r4. (XN’ + YN) - (XN + YN' ) 

^2W at P " 2XY 

T+ (XN' - YN' ) + (YN - XN) 


From the Figure A-3, XN' - YN' = XY 

and YN - XN » YZ 

where XZ is the perpendicular drawn on YN from X 


at P - 0 


XY -f YZ 
2XY 


This equation is valid for the right half of the tube, more- 
over upto 0 = a. 
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FIG. A -3 
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For 0 > = 0, till 0 = y in the other half. 

To find out YZ, the enlarged view of the triangle XYZ is 
considered, as shown in Figure A-4. 

From this figure, /.NPS^ = , /.YPO = 90° as OS’^ is 

normal to the surface at point P. 

XY is very small, so that YN and PN are almost para- 
llel, and as XZ is perpendicular to YN, so XZ is also perpendi- 
cular to PN, i.e. /.PX^X = 900 in the triangle PX^X and as OS^ 
is perpendicular to XY, ZXPS^ = 90o. 

/NPX = 900 -?^ 

/X^XP = 1800 _ [/NPX +ZPX^X1, = 


Now from the right angle triangle YZX, 

^ = sin i.e. YZ = XY sin W 


1 


yy 1 YZ ^ 

^2W at P "" XY^ 0 2XY = ^ 2X7 


= i (1 + sin ) 


.•. i (1 + sin for 0 < 0 < a 


(A.5) 


is found out considering Figure A-4, 


From triangle PS^N, 

From triangle PS 2 Nt 

(M = d + r2, MS 2 ^ 
OO 1 = r2 COS0, S 2 N 


= n - 

- - - (e 

- 2 

PS2 


= r 2 sin0, MN 
= MN - MSo ; 


+ 11/2 + 0 ) 

+ 0 ) 


= (W - L), 

(W - L) - r2 


(A.6) 


sin0 
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From expression of at P' , 
F 


2W at P' 


Lim sin $2 

X'Y' ^ 0 TT ^ 


^2W at P' ~ ~ ^ ^2^ 

From the same figure, /.OP'02 = (0 - u) 

M'P'N' = /OP'N’ - Z0P'02 

= . (13 ^ n) := 2Tt - (0 + 

From triangle N'M'P*, M'P'N' = | - ©2 
Equating these two, 

2ii - (0 + ^ 2 ) = I - ©2 


or, 5^2 = + ©2 

where, Sj = tan”^ ^ 


P'M' = M'02 + O 2 P' = d + r 2 + ^2 cos(0 - n) 

= d + r2 - ^2 COS0 

M'N' = MN' - MM' = L - r 2 sin(0 - it) = L + 

. d + r,5 - r^ COS0 

«2 = lA, slni8 — ^ 

Therefore all the expressions for F2 ^y( 0) different 
of 0 are listed below. 


(A.8) 


(A.9) 

r 2 sin0 

(A. 10) 
values 
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= ■5(1 + sin for 0 ^ 

= 0 , for a < 

= ■^(1 - sin y ^ 


where 




1 = I - (0^ + 0), 


and = tan 


-1 r ■*■ ^ 2 ^ ■ ^2 1 

(W - L) - r2 sin 0 


'P. = 1.511 + e 0 , 


and ©o = tan' 


-1 r ^ 2 ^ “ ^2 1 

L T . ^ ^ ~ ^ 


L + sin0 


0 < a 

0 < Y 

0 < 2it 
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APPENDIX B 

the LINE-BY-LINE METHOD 

A line-by-line method is a conventional combination 
of TDMA (Tri-Diagonal Matrix Algorithm) method for one dimen- 
sional situations and the Gauss-Seidel iterative method. Accor- 
ding to this method a grid line (in a particular direction) is 
chosen and the values of temperatures, T* s (or any dependent 
variable) along the neighbouring lines are assumed or known by’ 
their latest values, and the T's along the chosen line are 
solved by TDMA. This procedure is followed for all the lines 
in the same direction and then the technique is repeated till 
the required convergence is attained. Sometimes to make the 
convergence faster the same procedure is repeated for the 
lines in other direction. 

To visualize the line-by-line method more clearly the 
Figure B-1 is considered. The discretization equations for the 
grid points along a particular line (either in x-direction or 
in y-direction) are considered. These equations contain the 
temperatures at the grid points (shown by crosses) along the 
two neighbouring lines. If these temperatures are substituted 
from their latest values, the equations for the grid points 
(shown by dots) along the chosen line would look like one 
dimensional equations and could be solved by TDMA. This pro- 
cedure is carried out for all the lines in y-direction and 
may be followed by a similar treatment for x-direction. 
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The convergence of the line-by-line method is faster 
because the boundary condition information from the ends of 
the line is transmitted at once to the interior of the calcu- 
lation domain, no matter how many grid points lie along the 
line. The sweep direction (i.e. the sequence in which lines 
are chosen) plays an important role in the convergence of the 
iterative technique. A sweep from upstream to downstream 
would produce much faster convergence than a sweep against 
the stream. 
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APPENDIX C 

INPUT DATA TO THE PROGRAM 


i) Geometry of the tube: 


r^ (m) 

— r— 

i Th = 

t 

(r^ - r^ ), (m) 



0.005 

0.12 


0..01 

0.02 

0.03 

isition of the 

tube: 


d (m) 

1 W (m) 

1 

i L (m) 

t 

0.07 

0.4 

W 

2 

0.10 

1.0 

0.13 

20.0 


iii ) Conductivity of the tube material, ■): 

k = 35, 48 and 55. 

o 

iv) Known temperatures (®K) 

Tg = 420o0, Tg = 300.0, T^ = 830.0 

v) Radiation parameters: 

a = 0.567 X 10“"^ W/U^.k"^), e = 0.79 

vi ) Flow parameter: 

Re = 100, 1000, 1500 and 2000 

where, n(2r.) 

Re is defined as: Re = — 

vii) Number of grid points taken for numerical solution; 
Nj, = 50, N^,^ = 30, = 40 
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viii ) Properties of fluids! 


Pr . 

“T — 

i ■ 

I 

i w 

; m.K 

I 

t 

f 

1 

t 

j. 

a 

m2 

IB 

S 

X lo"^ 

“T 

f 

f V 

f 

! m2 

t 

f s 

i X 10^ 

1 

t 

51 

0.260 

0.932 

4.75 

84 

0.132 

0.663 

5.60 

175 

0.135 

0.710 

12.40 

490 

0.138 

0.769 

37.50 
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appendix D 


listing of the computer program 
NOMENCLATURE t 


A 

- coefficient of TDM 

! NR 


AK 


I NRF 

= Nrf 

AKF 


! NPR 

= Pr 

AL 

= L 

! PI 

= IT 

ALPHA 

« OC ■ 

! OB 

=5, 

ANU , 

== Nu 

I R1 

» r, 

ANUF 

= 0 

R2 


6 

« coefficient of TDM 

! REMOLD 

= Re 

C 

= coefficient of TDM 

! RI 

= rj 

D 

= coefficient of TDM 

\ SIG 

= <r 

D4 

= d 

1 T 

= T 

DPHI . 

= 

1 TFB 


DR 

= 41-5 

TFP 

= Tp 

DRF 

= 

! TS 

If 

DX 

= Ax 

! TW 

If 

—1 

2 

EPS 

= € 

! TUB 

= T, 

EPSN 

= Accuracy factor for 

! U 

= u 


convergence check 



F2W 


! UAVER 

rs U 

H 

= h 

! V 

= final 




vector 

1 

= i 

i W 

= U 


J . , = J 

K = k 


soln* 

ofTDMA 
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M* T E C H T H E S I S 

c title: conjugate heat transfer in a 

C LAMINAR P I P E FLOW SUBJECTED TO 

C non-uniform RADIATION FROM A LARGE 

C HEATED WALL : T H E CASE OF THERMALLY 

C U N D E U E L 0 P E D F L 0 W 

C main PROGRAM 

C *>l<^*)k)l<*;{<>r<)C>l«>l(*>«>K)k>K>K)K)K)K)l:***:(ok>K)K*>IC3K*)K*)|c>K)|c)|{)|{)K*)|cHc)lc)K>l<*3|c*)K>K)|c>tc>l<***>K>l<Hc>K 

C0MM0N/AREA1/D4 f R2 r AL y W !- F2W 
C0MM0N/AREA2/A y B y C y D » V 

DIMENSION T(50y50) yA<50) yB(50) y C ( 50 ) y D ( 50 ) y 0 ( 50 ) y F2W ( 50 ) 
DIMENSION TFPCSOySO) yTTl (50) yS(50) yU(50) yRI<50) yTWALL(50) 
DIMENSION Q(50) 

OPEN ( UNIT=--2:L y DEVICE-^ ' DSK " y FILE= ' MASTER . OUT ' ) 

AL-==0 * 5 y W= 1 ♦ 0 y R 1 =0 * 1 2 y R2==0 ♦ 1 4 y D4=0 ♦ 1 y 

AK=48 y 0 y TS==300 y 0 y SIG=0 * 567E-07 y TW=830 y 0 y EPS=0 y 79 

AKF=Oy 132yNRF==30 

NFM:l.-NRF-l 

ALPHA==0y663E-07 

RENOLD-1500y 0yNPR~84y ANUF=Oy056E-04 
UAMER==REN0LD*0 y 5>KANUF/R1 
TYPE *yUAOER 

WRITE ( 21 y 3333) NPRyRENOLDyANUFyUAMER 



3333 P'0I^'MAT(2Xj I6j2XyF6. 1 !-2Xf E8*3)-2Xi-F10t4/) 
DRF---R1/FL0AT(NRF-1) 

P 1 ==3 * 1 415927 i NPH I---=40 J NR==50 
TP I ==2 I 

dphi=:tp:i;/float ( nphi ) 

DR== ( R2"~R1 > /FLOAT (NR-NRF ) 

Pl==2 . :«Rl>KDRF+nRF)KDRF 
P2--=2 ♦ *Rl)HDR-DR>|cIiR 
R”Pl>lcAK/(P2)|«AKF) 

RR~< R>KDR*DR+nRF>l<i:iRF)/(DPHI*riPHI ) 

A1=:2*)KR1:(<R1 

A3~A1)KR 

A2=-<A3+Al-f2t,*RR) 

S2=2*>1«DR*SIG*EPS/AK 
S2:=S2* ( 2 ♦ *R2*R2+R2*DR ) 

DX-=1*0 

DRDF-DR/DPHI 

DRDPF=i:iRF/DPHI 

EPSN=-0 ♦ 1 

DO 10 I^lyNRF 

RI(I)=-DRF)K(I~1) 

10 U ( I ) "1 ♦ 0-RI ( I ) )KR.I ( I ) / ( R1*R1 ) 

DO 15 :C=-NRFflyNR 
1 5 ' R I ( I ) ==R 1 +DR* ( I “NRF ) 

DO 205 I-lyNRF 
DO 205 J=lyNPH;r 
205 TFP( I y J)“420*0 

CALL CONFIG(lyNPHI) 
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C0NST=::4 ^ 0*DF^F")KDRF- )i<UAVE:R/ ( ALPF^A*DX ) 

XX-=:0*0 

DO 202 K-=lx50 
DO 40 I-lxNR 
DO 40 J==lxNPHI 
40 T(I.xJ)=450* 

D099 ITER =1x200 
TOLD=-T(NRxl) 

C J=1 )|<)i<>|(***3K!|c)l<)K*)K)K)K)|c*)|c 

DO 1651 I=2xNRF'-l 
A ( I) =2 . *R I ( I ) *R I < I ) -R I ( I ) JkDRF 
C ( I ) =2 * JkRI ( I ) *RI ( I ) +RI ( I ) >KDRF 

B ( I ) =-4 * *RI< I) )KRI ( I ) “4 . *DRDPF*DRDPF-CONST*RI ( I ) JkRI ( I ) *U ( I ) 
1651 D ( I ) =-CONST)kRI ( I ))KRI < I ) *U ( I ) *TFP ( I x 1 ) -2 ♦♦DRDPFJKDRDPF* ( T ( I x 2 ) + 

IT(IxNPHI)) 

A( 1 )=0*0 

B ( 1 ) = 1 * 0 X C (1 ) = 1 ♦ 0 X D ( 1 ) = 0 ♦ 0 
D0101I=NRF+1 xNR 

A ( I ) =2 * >I<RI ( I ) 5KRI ( I )""RI ( I )*DR , 

C (I ) =2 ♦ skRI < I ) JKRI < I ) +DR>KRI ( I ) 

B ( I ) =-4 . )|CRI ( I ) >KRI < I ) “4 * *DRDP>kDRDP 
1 01 D ( I ) =-2 X )(<DRDP)KDRDP*T (I x NPHI ) -2 * *DRDP)kDRDP*T (1x2) 

A(NR)=4*>kR25l<R2 
C(NR)=0*0 

B ( NR ) =B ( NR ) -"4 X *82*1 ( NR X 1 ) *5k3 X 

D ( NR ) =D< NR ) -•S2*F2U ( 1 ) *TU**4x-S2*( 1 x -F2W ( 1 ) )*TS**4x -3x *82 
l*T(NRxl)**4, 

A(NRF)=AlxB(NRF)=A2xC(NRF)=A3 
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102 

C 


210 


104 


D ( NRF" ) =-"-F^R)KT < NRF y 2 ) ~RR>kT < NRF y NPFII ) 

CALL TRiriAGdyNR) 

DO 102 1=1 y NR 

T(;Lyl)=0(.T) 

>K)K)K)K)|()K*Hc)|()K)|C)|c)K>|c>(:)(:)|c)K>|c>k)((>K>K J=2yNPFlI“l )|c>k)|c)i<5K>|c*>K)K>K>|c*!K>K5K5K)K3K)|c)K)K5K>K* 

DO 103 J=2yNPHI-l 
DO 210 1=2 y NRF- 1 

A ( I ) ■•■“2 ♦ )kR I ( I ) *R I ( I) -R Id) )KDRF 

C d ) =2 ♦ )KR I (I ) )KR I ( I ) +RI d ) >»;DRF 

B ( I ) =-4 . *RI d ) )KRI d ) -4 ♦ *DRDPF>KDRDPF-CONST*RI d ) *RI< I ) *U d ) 

D (I ) =-CONST)KRI d ) *RI d ) *0 d ) >l!TFP d y J ) -2 « )|{DRDPF>KDRDPF>lc ( T d y J+1 ) 4 
ITdy J-1) ) 

Ad)=0.0yC<l):=l40 
Bd)=-l,0yDd)==0*0 
DO 104 I=NRF+lyNR 
A d ) =2 . *R I d ) :^RI d ) “RI d ) JkDR 
C d ) ===2 ♦ )KRI d ) JKRI d ) 4-RI d ) !lcDR 
BCD =-4 * )KRI d ) *RI d ) -4 * )ki;iRDP*DRDP 

D d ) =-2 ♦ )KDRDP>KDRDP>kT d y J-1 ) -2 * JKDRDPJKDRDP*! < I y J4-1 ) 

A(NR)^=4*5KR2*R2 

C(NR)=0»0 

B(NR)=B(NR)~4y*S2)KT(NRy J)**3. ’ 

D ( NR ) =D < NR ) -S2>KF2W ( J ) )KTW)ic5K4 » -S2* C 1 * ~F2W ( J ) ) tTSttA . -3 ♦ *S2 
l)KT(NRy J)*)K4y 

A ( NRF ) - A 1 y B ( NRF ) ~A2 y C ( NRF ) 

D CNRF)=-RR5t: CT (NRFy J+l )4-T(NRFy J-1 ) ) 

CALL TRIDAGdyNR) 


DO 105 1=1 y NR 
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105 T(I»J)-V(I) 

103 CONTINUE 

DO 215 I=2yNRF-l 
A ( I ) -2 ♦ *R I ( .1 ) *R I ( I ) -R I ( I ) !j<DRF 
C ( I ) -2 . >I<R;C < I ) )|{RI < .T ) -f-RI ( I ) !{{DRF 

B ( I ) =-4 . *R.T( I ) tR 1 ( I ) ~4 4 *DRDPF*DRriPF“CONST*RI ( I) )tcRI ( I ) )lcU ( I ) 

21 5 D ( I ) =- CONST *RI ( I ) tRl ( I ) 5|cU ( I ) *TFP ( I y NPHI ) -2 . *DRDPF*DRDPF* ( T ( I y 1 ) 

l+T(IyNPHI”-l) ) 

A ( 1) =0 4 0 y B ( 1 ) “:-l 4 0 y C ( 1 ) =1 4 0 y D ( 1 ) =0 4 0 
no 106 I-NRF+lyNR 
A ( I) “2 4 )KR I ( I ) )KRI < I ) “R I ( I ) *DR 
C ( I ) =--2 4 )KRI < I ) *RT. ( I ) i RI ( I )*DR 
B ( I ) = "4 4 *RI ( I ) *RI ( I ) -4 4 )|cnRnP>lcriRDP 

106 n ( I ) ==“2 4 *nRDP*nRBP)KT ( I y NPHI-1 ) -2 4 >KnRDP*DRriP*T ( I y 1 ) 

A<NR)=4 4*R2>KR2 

C(NR)=040 

B ( NR ) =B ( NR ) -4 4 *S2*T ( NR y NPHI ) **3 4 

n < NR ) =D ( NR ) -S2)KF2W < NPHI ) *TU*3K4 » -S2* ( 1 * -F2W ( NPHI ) ) *13**4 4 -32*3 4 * 
IT ( NR y NPHI >**44 
A ( NRF ) =A 1 y B < NRF ) ==A2 y C ( NRF > =A3 
n (NRF ) :--RR*T ( NRF y 1 ) -RR*T ( NRF y NPHI-1 ) 

CALL. TRIDAG(lyNR) 
no 107 1=1 y NR 

107 T(IyNPHI>=U(I) 

IF(AB3(TOLn-T<NRyl)>*LE4EPSN>GO TO 110 

99 CONTINUE 


WRITE(21y 115) 



115 FORMAT (2X» 'NOT CONVEF{GEDV/) 

WRITE (21 1-777) 

777 F0RMAT(2Xy 'INTER FACE TEMF-v y J=--l y NF'HI ' / ) 

WRITE(21y779) (TCNRFy J) yJ=lyNPHI) 

779 FORMAT (8(F12* 3) ) 

WRITE ( 21 y 885) 

885 FORMAT (2Xy 'TEMP* FOR THE WHOLE REGION WITH J^lyNPHI 

lJRITE(21y!K) ( (TCIyJ) yJ==lyNPHIy5) yI = lyNR) 
no 889 I==lyNRF 
no 889 J=lyNPHI 

889 TFP(IyJ)===T(IyJ) 

C H E A T T R A N S F E R C A L C U L A T IONS 

no 11 I=lyNRF 
no 12 J=lyNPHI 

12 TTl ( J)==T<Iy J) 

TT1(NPHI + 1)==TT1(1) 

CALL SIMPSN<OyTPIyNPHIyTTlyXINTGL) 
S(I)==XINTGL>KRI(I)#U(I) 

11 CONTINUE 

171 CALL SIMPSN(OyRlyNFMlySyXINTGL) 

TFB=2 * 3|cX I NTGL/ ( PI >KR 1 *R 1 ) 
no 13 J-1 tNPHI 
TWALL( J)~T<NRFy J) 

13 Q ( J ) - ( T ( NRF-f 1 y J ) “T ( NRF y J ) ) /HR 
TWALL ( NPHIfl ) “TWALL ( 1 ) 

Q(NPHIfl )==Q( 1 ) 

CALL SIMPSN(OyTPIyNPHIyTWALLyXIN‘TGL) 


TWB"XINTGL/TPI 
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CALL SIMF'SN(0»TPI»NPHI i-Qi-XXNTGL) 

aB==X.tNTGL/TPI 

H“AKH<QB/<TU)B-TFB) 

ANU=2*)KH>KFLI./AKF 
TYPE *»ANU 
XX--=XX+DX 

WRI TE ( 21 y 1 4 ) XX V TUB 1 - QB y TFB H"l 1 - ANU 

14 FORMAT < 2X y F8 ♦ 2 y 4X y F8 ♦ 3 y 4X y F7 ♦ 2 y 4X y F8 ♦ 3 y 4X y F7 ♦ 3 y 4X y F8 ♦ 3 ) 

202 CONTINUE 

STOP 
END 

C S U B R 0 U T I N E : T R I D A G 

C (THIS SUBROUTINE SOLOES SIMULTANEOUS LINEAR NON-HOMOGENIOUS EONS* 

C BY TRI-DIAGONAL MATRIX SOLUTION) 

SUBROUTINE TRIDAG(IFyL) 

C0MM0N/AREA2/AyByCyDy V 
D I MENS ION BETA (S 1 ) y GAMMA (51) 

DIMENSION A(50) yB(50) yC(50) yD(50) y0(50) 

C COMPUTE INTERMEDIATE ARRAYS BETA AND GAMMA 

BETA( IF)=:B<IF) 

GAMMA(IF)==D(IF)/BETA(IF) 

IFPl^IFTl 
DO 1 I=0;FPlyL 

BETA( I )=~B( D-AC I )*C( I"1 )/BETA( I-l ) 
r gamma ( I >- (D( I )-A( I )>}:GAMMA ( I~1 ) )/BETA( I ) 
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C COMPUTE FINAL SOLUTION VECTOR V 

V (L) “GAMMA (L) 

LAST==L-IF 
DO 2 K=li-LAST 
I=L-K 

2 V(I)=GAMMA(I)--C(I)>KV(I + :L)/BETA(I) 

RETURN 

END 

C S U B R 0 U T I N E : C 0 N F I G 

C SUBROUTINE TO CALCULATE CONFIGURATION FACTOR 

jjSi /fv /jp. /fw wY^ ^Y** ^ NT* m' ^Y ^Y ^ ^Y ^Y ‘Y ^Y •Y ‘Y ^Y Y* ‘Y "Y Y^ 'Y Y^ Y^ Y* Y^ Y* Y^ ^Y ^Y 

SUBROUTINE CONFIG ( IS y NPHI) 

COMMON/AREAI /D4 y R2 y AL y W y F2W 
DIMENSION F2UK50) 

PI==3* 1415927 
TPI=2y)|{PI 

Dl=SQRT(D45KD4+2y3|cR2)KD4) 

BETA=ATAN<D1/R2) 

Bl=2*>l<R2)k<D4+R2) 

C1=R2*R2“(W~AL>*(U)"-AL) 

A1“(D4+R2)>K(D4+R2)K W-AL)*(U“AL) 

D1==SQRT(B1>KB1“4**A1*C1 ) 

D2==4*>KAl*Al--(Dl+Bl)>K(DlfBl) 

D2==SQRT(D2) 

ALPHA1=ATAN(D2/(B1+D1) ) 
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D2~-4*)KA1*A1-"(D1~B1 )>K(ril"-Bl ) 

D2:==SQRT(Li2) 

ALPHA2=ATAN < Ii2/ ( Bl-Dl ) ) 

IF ( ALPHA2 . LT ♦ 0 ♦ 0 ) ALPHA2~ALPHA2+PI 
IF(ALPHA2.LT»Bli::TA) GO T05 
ALPHA=-Al..PHA2 
GO TO 6 

5 ALPHA:=ALPHA1 

A C1==R2)I<R2-AL)KAL 

A1=AL:KAL+ ( ri4+R2 ) * ( D4+R2 ) 

D 1 ^-SQR T ( B 1 *B 1 ”4 . * A;l. *C 1 ) 
D2-4**Al>KAl-(Dl+Bl>)K<i:il+Bl) 

ri2-=SQRT(D2) 

/ 

ALPHAl=ATAN(D2/<ri:l fBl) ) 

B 2 = 4 . * A 1 * A 1 ~ ( D 1 - B 1 ) * ( D 1 •“ B 1 ) 

D2==SQRT(D2) 

ALPHA2=-=ATAN(D2/(B1"D1) ) 

IF ( ALPHA2 ♦ LT « 0 ♦ 0 ) ALPHA2==ALPHA2+PI 
IF(ALPHA2J...T*BETA) GO TO 7 
DELTA--=ALPHA2 
GO TO 16 

7 DELTA-=ALPHA1 

16 nPHI=TPI/FLOAT (NPHI ) 

6AMA=TPI -DELTA 

PHI=0«0 

K--=l 

10 X=(D4+R2-R25KC0S(PHI))/(W-AL~R2*SIN(PHI)) 

THETA=ATAN(X) 


I 
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SYl:=0.5)kPI- THETA-PHI 
F1W-0»5*(1«+SIN(SY1) ) 

F2W(K)=---F1W 

WFi;ITE(21y*)PHIyFlUI 

PHI==PHI+DPH1 

K=K+1 

IF ( PH I, LE» ALPHA) GO TO 10 

15 riw==o*o 

WFMTE(21i-*)PHI»FlW 

F2W(K)-F1W 

K=::Kfl 

PHI=--PHI+i;iPHI 
IF ( PHI ♦LE. GAMA) GO TO 15 
20 X=(D4+R2-r-;:2>KC0S(PHI) )/<AL+R2*SIN<PHI)) 

THETA=ATAN<X) 

SY1-1.5>KPI + THETA-PHI 
F1U=0*5*(1«-SIN(SY1)) 

WRITE(21f)ic)PHIi-FlW 

F2W<K)=:F1W 

K=K+1 

PHI=PHI+DPHI 
IF(PHI,LE*TPI)GO TO 20 
RETURN 
END 
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THIS SUBPROGRAM INTEGRATES ANY GIVEN FUNCTION OVER A FINITE 
LIMIT USING SIMPSON'S I/S-RD RULE 
jK )K XOK ** ** )!<)K>lc >K )K*>K ****** )K *** He >K *** >K * ^ ****** **>!<* *>1! >lc ***** )lc:i<*5|c* * 

SUBROUTINE SIMPSN ( XMIN y XMAX » Nr DUMMYF > XINTGL ) 

DIMENSION DUMMYFOO) 

H= ( XMAX-XM IN ) /FLOAT ( N ) 

SUM==0.0 
DO 4 1=2 rN 
IF(M0D<Iy2) )2y2y3 
SUM=SUM+4,*DUMMYF< I ) 

GO TO 4 

SUM=SUM+2 ♦ *DUMMYF ( I) 

CONTINUE 

XINTGL=H/3 ♦ * ( DUMMYF ( 1 ) +SUM+DUMMYF ( N+1 ) ) 

RETURN 

END 


******************** E N D **************************** 



